{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Datashading LandSat8 raster satellite imagery\n",
    "\n",
    "Datashader is fundamentally a rasterizing library, turning data into rasters (image-like arrays), but it is also useful for already-rasterized data like satellite imagery.  For raster data, datashader uses the separate [xarray](http://xarray.pydata.org/) library to re-render the data to whatever new bounding box and resolution the user requests, and the rest of the datashader pipeline can then be used to visualize and analyze the data.  This demo shows how to work with a set of raster satellite data, generating images as needed and overlaying them on geographic coordinates using [HoloViews](http://holoviews.org) (version >1.9), [GeoViews](http://geo.holoviews.org), and [Bokeh](http://bokeh.pydata.org)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "import geoviews as gv\n",
    "import datashader as ds\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "from holoviews.operation.datashader import regrid, shade\n",
    "from bokeh.tile_providers import STAMEN_TONER\n",
    "\n",
    "hv.extension('bokeh', width=80)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load LandSat Data \n",
    "\n",
    "LandSat data is measured in different frequency bands, revealing different types of information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "band_info = pd.DataFrame([\n",
    "        (1,  \"Aerosol\", \" 0.43 - 0.45\",    0.440,  \"30\",         \"Coastal aerosol\"),\n",
    "        (2,  \"Blue\",    \" 0.45 - 0.51\",    0.480,  \"30\",         \"Blue\"),\n",
    "        (3,  \"Green\",   \" 0.53 - 0.59\",    0.560,  \"30\",         \"Green\"),\n",
    "        (4,  \"Red\",     \" 0.64 - 0.67\",    0.655,  \"30\",         \"Red\"),\n",
    "        (5,  \"NIR\",     \" 0.85 - 0.88\",    0.865,  \"30\",         \"Near Infrared (NIR)\"),\n",
    "        (6,  \"SWIR1\",   \" 1.57 - 1.65\",    1.610,  \"30\",         \"Shortwave Infrared (SWIR) 1\"),\n",
    "        (7,  \"SWIR2\",   \" 2.11 - 2.29\",    2.200,  \"30\",         \"Shortwave Infrared (SWIR) 2\"),\n",
    "        (8,  \"Panc\",    \" 0.50 - 0.68\",    0.590,  \"15\",         \"Panchromatic\"),\n",
    "        (9,  \"Cirrus\",  \" 1.36 - 1.38\",    1.370,  \"30\",         \"Cirrus\"),\n",
    "        (10, \"TIRS1\",   \"10.60 - 11.19\",   10.895, \"100 * (30)\", \"Thermal Infrared (TIRS) 1\"),\n",
    "        (11, \"TIRS2\",   \"11.50 - 12.51\",   12.005, \"100 * (30)\", \"Thermal Infrared (TIRS) 2\")],\n",
    "    columns=['Band', 'Name', 'Wavelength Range (µm)', 'Nominal Wavelength (µm)', 'Resolution (m)', 'Description']).set_index([\"Band\"])\n",
    "band_info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_path = '../data/MERCATOR_LC80210392016114LGN00_B%s.TIF'\n",
    "bands = list(range(1, 12)) + ['QA']\n",
    "bands = [xr.open_rasterio(file_path%band).load()[0] for band in bands]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rendering LandSat data as images\n",
    "\n",
    "The bands measured by LandSat include wavelengths covering the visible spectrum, but also other ranges, and so it's possible to visualize this data in many different ways, in both true color (using the visible spectrum directly) or false color (usually showing other bands).  Some examples are shown in the sections below.\n",
    "\n",
    "### Just the Blue Band\n",
    "\n",
    "Using datashader's default histogram-equalized colormapping, the full range of data is visible in the plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "opts.defaults(opts.RGB(width=600, height=600))\n",
    "\n",
    "nodata= 1\n",
    "\n",
    "def one_band(b):\n",
    "    xs, ys = b['x'], b['y']\n",
    "    b = ds.utils.orient_array(b)\n",
    "    a = (np.where(np.logical_or(np.isnan(b),b<=nodata),0,255)).astype(np.uint8)    \n",
    "    col, rows = b.shape\n",
    "    return hv.RGB((xs, ys[::-1], b, b, b, a), vdims=list('RGBA'))\n",
    "\n",
    "tiles = gv.WMTS(STAMEN_TONER)\n",
    "tiles * shade(regrid(one_band(bands[1])), cmap=['black', 'white']).redim(x='Longitude', y='Latitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will usually want to zoom in, which will re-rasterize the image if you are in a live notebook, and then re-equalize the colormap to show all the detail available.  If you are on a static copy of the notebook, only the original resolution at which the image was rendered will be available, but zooming will still update the map tiles to whatever resolution is requested.\n",
    "\n",
    "The plots below use a different type of colormap processing, implemented as a custom transfer function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.utils import ngjit\n",
    "nodata= 1\n",
    "\n",
    "@ngjit\n",
    "def normalize_data(agg):\n",
    "    out = np.zeros_like(agg)\n",
    "    min_val = 0\n",
    "    max_val = 2**16 - 1\n",
    "    range_val = max_val - min_val\n",
    "    col, rows = agg.shape\n",
    "    c = 40\n",
    "    th = .125\n",
    "    for x in range(col):\n",
    "        for y in range(rows):\n",
    "            val = agg[x, y]\n",
    "            norm = (val - min_val) / range_val\n",
    "            norm = 1 / (1 + np.exp(c * (th - norm))) # bonus\n",
    "            out[x, y] = norm * 255.0\n",
    "    return out\n",
    "\n",
    "def combine_bands(r, g, b):\n",
    "    xs, ys = r['x'], r['y']\n",
    "    r, g, b = [ds.utils.orient_array(img) for img in (r, g, b)]\n",
    "    a = (np.where(np.logical_or(np.isnan(r),r<=nodata),0,255)).astype(np.uint8)    \n",
    "    r = (normalize_data(r)).astype(np.uint8)\n",
    "    g = (normalize_data(g)).astype(np.uint8)\n",
    "    b = (normalize_data(b)).astype(np.uint8)\n",
    "    return hv.RGB((xs, ys[::-1], r, g, b, a), vdims=list('RGBA'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### True Color\n",
    "\n",
    "Mapping the Red, Green, and Blue bands to the R, G, and B channels of an image reconstructs the image as it would appear to an ordinary camera from that viewpoint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_color = combine_bands(bands[3], bands[2], bands[1]).relabel(\"True Color (R=Red, G=Green, B=Blue)\")\n",
    "tiles * regrid(true_color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the raster data will only refresh to a new resolution if you are running a live notebook, because that data is not actually present in the web page; it's held in a separate Python server.\n",
    " \n",
    "### False Color\n",
    "\n",
    "[Other combinations](https://blogs.esri.com/esri/arcgis/2013/07/24/band-combinations-for-landsat-8/) highlight particular features of interest based on the different spectral properties of reflectances from various objects and surfaces, with full data redrawing on zooming if you have a live Python process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combos = pd.DataFrame([\n",
    "    (4,3,2,\"True color\",\"\"),\n",
    "    (7,6,4,\"Urban\",\"False color\"),\n",
    "    (5,4,3,\"Vegetation\",\"Color Infrared\"),\n",
    "    (6,5,2,\"Agriculture\",\"\"),\n",
    "    (7,6,5,\"Penetration\",\"Atmospheric Penetration\"),\n",
    "    (5,6,2,\"Healthy Vegetation\",\"\"),\n",
    "    (5,6,4,\"Land vs. Water\",\"\"),\n",
    "    (7,5,3,\"Atmosphere Removal\",\"Natural With Atmospheric Removal\"),\n",
    "    (7,5,4,\"Shortwave Infrared\",\"\"),\n",
    "    (6,5,4,\"Vegetation Analysis\",\"\")],\n",
    "    columns=['R', 'G', 'B', 'Name', 'Description']).set_index([\"Name\"])\n",
    "combos\n",
    "\n",
    "def combo(name):\n",
    "    c=combos.loc[name]\n",
    "    return regrid(combine_bands(bands[c.R-1],bands[c.G-1],bands[c.B-1])).relabel(name)\n",
    "\n",
    "layout = combo(\"Urban\") + combo(\"Vegetation\") + combo(\"Agriculture\") + combo(\"Land vs. Water\")\n",
    "\n",
    "layout.opts(\n",
    "    opts.RGB(width=350, height=350, xaxis=None, yaxis=None, framewise=True)).cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the various ways of combining aggregates supported by [xarray](http://xarray.pydata.org) are available for these channels, making it simple to make your own custom visualizations highlighting any combination of bands that reveal something of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Revealing the spectrum\n",
    "\n",
    "The above plots all map some of the measured data into the R,G,B channels of an image, showing all spatial locations but only a restricted set of wavelengths. Alternatively, you could sample across all the measured wavelength bands to show the full spectrum at any given location:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band_map = hv.HoloMap({i: hv.Image(band) for i, band in enumerate(bands)})\n",
    "\n",
    "def spectrum(x, y):\n",
    "    try: \n",
    "        spectrum_vals = band_map.sample(x=x, y=y)['z'][:-1]\n",
    "        point = gv.Points([(x, y)], crs=ccrs.GOOGLE_MERCATOR)\n",
    "        point = gv.operation.project_points(point, projection=ccrs.PlateCarree())\n",
    "        label = label = 'Lon: %.3f, Lat: %.3f' % tuple(point.array()[0])\n",
    "    except:\n",
    "        spectrum_vals = np.zeros(11)\n",
    "        label = 'Lon: -, Lat: -'\n",
    "    \n",
    "    return hv.Curve((band_info['Nominal Wavelength (µm)'].values, spectrum_vals), label=label,\n",
    "                    kdims=['Wavelength (µm)'], vdims=['Luminance']).sort()\n",
    "\n",
    "# x and y give the location in Web Mercator coordinates\n",
    "spectrum(x=-9880000, y=3570000).opts(width=800, height=300, logx=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now combine these two approaches to let you explore the full hyperspectral information at any location in the true-color image, updating the curve whenever you hover over an area of the image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tap = hv.streams.PointerXY(source=true_color)\n",
    "spectrum_curve = hv.DynamicMap(spectrum, streams=[tap]).redim.range(Luminance=(0, 30000))\n",
    "\n",
    "layout = tiles * regrid(true_color) + spectrum_curve\n",
    "layout.opts(\n",
    "    opts.Curve(width=450, height=450, logx=True),\n",
    "    opts.RGB(width=450, height=450))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Of course, just as for the raster data resolution, the plot on the right will update only in a live notebook session, because it needs to run Python code for each mouse pointer position.)\n",
    "\n",
    "As you can see, even though datashader is not a GIS system, it can be a flexible, high-performance way to explore GIS data when combined with HoloViews, GeoViews, and Bokeh."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
